% plot individual fit parameters for all curves

tunneling = 2 * pi * 925; %Hz
% compare these to the full parameters
full_data_file = 'expt_hydro_fit_summary.txt';
full_data = dlmread(full_data_file, '\t', 1, 0);
temps_full = full_data(:, 1);
[temps_full, index] = sort(temps_full);
temp_uncs_full = full_data(index, 2);
gammas_full = full_data(index, 7) * 1e6/tunneling;
gamma_errs_full = full_data(index, 8) * 1e6/tunneling;
ds_full = full_data(index, 9) * 1e6/tunneling;
d_errs_full = full_data(index, 10) * 1e6/tunneling;

files = dir('*/*/*_struct.mat');

% period_bin_endpts = [0,7,11,17,25]; %[0, 7]; %, 14, 17, inf];
period_bin_endpts = [0, 7, 14];
nbins = length(period_bin_endpts) - 1;

ds_stack = [];
for ii = 1:length(files)
    a = load(fullfile(files(ii).folder, files(ii).name));
    dsTemp = DiffuseSet();
    dsTemp.loadStruct(a.AsStruct, '', 0);
    ds_stack = cat(1,ds_stack,dsTemp);
end

periods = [ds_stack.Period];
period_uncs = [ds_stack.PeriodUnc];
temps = [ds_stack.Temp];
temp_uncs = [ds_stack.TempUnc];

fit_params = zeros(length(ds_stack), 3);
std_errs = zeros(length(ds_stack), 3);
for ii = 1:length(ds_stack)
    fit_params(ii,:) = ds_stack(ii).Fpdo(1:3);
    std_errs(ii,:) = ds_stack(ii).SEdo(1:3);
end
% fit_params = vertcat(ds.Fpdo);
% std_errs = vertcat(ds.SEdo);

gammas = fit_params(:, 1) * 1e6/tunneling;
gamma_errs = std_errs(:, 1) * 1e6/tunneling;

ds = fit_params(:, 2) * 1e6/tunneling;
d_errs = std_errs(:, 2) * 1e6/tunneling; 

plot_handles = [];
leg = {};

figure;
ax1 = subplot(1, 2, 1);
ax2 = subplot(1, 2, 2);
for ii = 1:nbins
    lowp = period_bin_endpts(ii);
    highp = period_bin_endpts(ii + 1);
    
    temps_bin = temps(periods < highp & periods > lowp);
    temp_uncs_bin = temp_uncs(periods < highp & periods > lowp);
    gammas_bin = gammas(periods < highp & periods > lowp);
    gamma_errs_bin = gamma_errs(periods < highp & periods > lowp);
    ds_bin = ds(periods < highp & periods > lowp);
    d_errs_bin = d_errs(periods < highp & periods > lowp);
    
	if ~isempty(temps_bin)

        ph = errorbar(ax1, temps_bin, gammas_bin, gamma_errs_bin, gamma_errs_bin, temp_uncs_bin, temp_uncs_bin, 'o');
        errorbar(ax2, temps_bin, ds_bin, d_errs_bin, d_errs_bin, temp_uncs_bin, temp_uncs_bin, 'o', 'Color', ph.Color);
        leg{end+1} = sprintf('Periods = %0.1f - %0.2f', lowp, highp);
        plot_handles(end+1) = ph;

        hold(ax1, 'on');
        hold(ax2, 'on');
    end
        
end

% finally, plot full data
ph = errorbar(ax1, temps_full, gammas_full, gamma_errs_full, gamma_errs_full, temp_uncs_full, temp_uncs_full, 'o');
errorbar(ax2, temps_full, ds_full, d_errs_full, d_errs_full, temp_uncs_full, temp_uncs_full, 'o', 'Color', ph.Color);

leg{end + 1} = 'Simultaneous Fit';
plot_handles(end + 1) = ph;

grid(ax1);
grid(ax2);
xlabel(ax1, 'T (t)');
xlabel(ax2, 'T (t)');
ylabel(ax1, 'Gamma (t/hbar)');
ylabel(ax2, 'D (a^2 * t/hbar)');
% ylim(ax1, [0, 2]);
ylim(ax2, [0, 6]);

legend(plot_handles, leg);